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Motivated by the observation of highly unstable flowing states in suspensions of microtubules and 
kinesin, we analyze a model of mutually-propelled filaments suspended in a solvent. The system 
undergoes a mean-field isotropic-nematic transition for large enough filament concentrations when 
the nematic order parameter is allowed to vary in space and time. We analyze the model in two 
contexts: a quasi-one-dimensional channel with no-slip walls and a two-dimensional box with pe- 
riodic boundaries. Using stability analysis and numerical calculations we show that the interplay 
between non-uniform nematic order, activity, and flow results in a variety of complex scenarios that 
include spontaneous banded laminar flow, relaxation oscillations, and chaos. 



I. INTRODUCTION 

Active hydrodynamics describes the collective motion of 
microscopic particles constantly maintained out of equi- 
librium by internal energy sources. Colonies of swarming 
bacteria, in vitro mixtures of cytoskeletal filaments and 
motor proteins and vibrated granular rods are common 
examples of active systems and now active has become 
standard terminology for any system whose constituents 
drive themselves mechanically by extracting and dissi- 
pating energy from their environment. Originating from 
pioneering works by Pedlcy and Kessler [1], Vicsek et 
al. [2] and Simha and Ramaswamy [3], active matter re- 
search has blossomed to encompass diverse systems and 
scales ranging from animal groups to subcellular matter 
[4]- 

Because active particles typically have elongated 
shapes, their collective behavior has been often described 
using the language of liquid crystals [5, 6]. In this re- 
gard, an important distinction among active particles 
concerns the possibility of forming phases characterized 
by nematic or polar order. While all elongated parti- 
cles can form a nematic phase at sufficient densities and 
levels of activity, particles which have an asymmetry as- 
sociated with their mutual interaction can additionally 
form a phase characterized by a non-zero macroscopic 
polarization. Active particles can also be distinguished in 
terms of their locomotion characteristics: self-propelled 
particles (SPP) are endowed with an internal engine and, 
typically but not necessarily, with appendages that al- 
low them to swim in a fluid or crawl on a substrate. 
For example, bacteria [7], large animals such as fish or 
birds [8], and catalytic motors [9] belong in this cate- 
gory. Cytoskeletal filaments, on the other hand, cannot 
propel themselves, but move in a solvent through the 
action of motor proteins, which are themselves powered 
by the hydrolysis of adenosine triphosphate (ATP). Bun- 



* Electronic address: lgiomi@seas.harvard.edu 
t Electronic address: hagan@brandeis.edu 



dies of molecular motors attach to pairs of filaments and, 
during an ATP cycle, slide the filaments with respect 
to each other. We will refer to this type of active ele- 
ments as mutually-propelled particles (MPP). There is 
finally a third class of systems in which activity is pro- 
vided through vibration. Vertically shaken granular rods, 
for instance, gain and dissipate energy while bouncing on 
a substrate, resulting in a two-dimensional motion along 
their major axis [10, 11]. 

Most theoretical effort in modeling active systems 
characterized by liquid crystalline order has focused on 
constructing hydrodynamic equations that, in addition 
to the usual liquid crystalline elasticity, can account for 
the additional forces and currents originated from the 
activity. This task has been achieved by incorporating 
phenomenological non-equilibrium terms in the hydro- 
dynamic equations of nematic and polar liquid crystals 
[3, 6, 12] or by applying the tools of non-equilibrium 
statistical mechanics to specific microscopic models [13- 
16]. This program has generated a variety of predictions, 
which include the existence of giant density fluctuations 
in active nematics [10, 17, 18], spontaneously flowing 
states [12, 19-25], unconventional rheological properties 
[24, 26-30] and a plethora of a novel hydrodynamic in- 
stabilities with no counterpart in passive complex fluids 
[3, 16, 31-34]; a recent overview can be found in [4]. 

In spite of the vast theoretical work on active liquid 
crystals, little consideration has been given to the pos- 
sibility of spatial and temporal variations in the order 
parameter; a recent exception is the work of Mishra et 
al. [35] who considered self-propelled polar rods mov- 
ing on a frictional substrate, with a density-driven mean 
field transition from the isotropic to the polar phase and 
showed that when the self-propulsion velocity exceeds a 
threshold value, the uniformly polarized moving state be- 
comes unstable to spatial fluctuations which organize into 
stripes of different density and polarization. 

Here we consider the case of an active nematic suspen- 
sion motivated by the observation of highly unstable flow- 
ing states in assemblies of microtubules and kinesin [36], 
a model for mutually-propelled elongated particles in a 
solvent. The system undergoes a mean-field isotropic- 



2 



nematic transition for large enough filament concentra- 
tions and the nematic order parameter is allowed to vary 
in space and time. We use stability analysis and numer- 
ical simulations to analyze the model in two geometries: 
a quasi-one-dimensional channel with no-slip walls and 
a two-dimensional box with periodic boundaries. In the 
channel geometry, moderate activity levels lead to spon- 
taneous laminar flow, as seen in earlier works [12, 19] 
that assumed a constant magnitude of the nematic order 
parameter. Upon increasing the activity past a threshold 
value, however, fluctuations in magnitude of the nematic 
order parameter leads to oscillatory flow in which the ne- 
matic director periodically switches orientation. In the 
two-dimensional box, the interplay between non-uniform 
nematic order, activity, and flow results in a variety 
of complex scenarios that include spontaneous laminar 
flow, relaxation oscillations reminiscent of excitable me- 
dia, and chaos. A detailed analysis allows us to uncover 
the origin of oscillations in the system and characterize 
the chaotic regime, wherein we see behavior consistent 
with turbulent flow even in the low Reynolds number 
regime, expanding on and complementing a recent short 
report of some of our findings [37]. 

This article is organized as it follows. In Sec. II we 
introduce the hydrodynamic equations for an active sus- 
pension of mutually propelled filaments. In Sec. Ill we 
analyze the equations for a quasi-one-dimensional chan- 
nel of infinite length and finite width endowed with no- 
slip walls. In Sec. IV we consider an active nematic 
suspension in a two-dimensional container with periodic 
boundaries. We then present a minimal model that 
demonstrates oscillatory behavior, and characterize the 
the chaotic regime. Finally, we present our conclusions 
in Sec. V. 



II. HYDRO DYNAMICAL EQUATIONS OF 
MOTION 

A fluid of orientable fore-aft symmetric particles can 
generally exist in two phases: isotropic (I) and nematic 
(N). In the latter phase, the particles are orientationally 
ordered with an average orientation characterized by the 
nematic director field n. For microscopic particles in sus- 
pension, such as colloidal rods or biological filaments, the 
IN transition is driven by density: when the concentra- 
tion of particles overcomes some critical value c*, the 
particles form a nematic phase in order to maximize en- 
tropy. In a two-dimensional equilibrium fluid of slender 
rods the critical concentration is given by c* = 3ir/2£ 2 
where £ is the length of the rods [38] and the phase tran- 
sition is of the Kosterlitz-Thouless type. The anisotropy 
of a nematic phase is expressed through the nematic ten- 
sor Qij [39], which for uniaxial nematics reads: 

Qij = S (niUj - i Sij) (1) 



The nematic tensor Qij is by construction traceless and 
symmetric, thus in d = 2 it consists of only two indepen- 
dent degrees of freedom. The nematic phase has orienta- 
tional order (S ^ 0) and is invariant under inversion of 
the director field: n — > — n. Here the extent of nematic 
alignment is expressed in terms of a scalar nematic order- 
parameter S: 

<?=^l_<d|a.n| 2 -l), (2) 

where a is the axis of the molecules, d is the dimension 
of the system and the angular brackets denote a thermal 
average. In a suspension of rod-like particles, S depends 
on the local concentration of the particles and, in equilib- 
rium passive systems, is constant across the sample since 
diffusion drives the fluid toward a homogeneous state. In 
active systems, however, activity can build up density in- 
homogeneities and the order parameter may exhibit spa- 
tial fluctuations. Moreover, since the effects of activity 
are generally enhanced by local orientational order, cou- 
pling between order, activity and flow can amplify these 
fluctuations. In the following we describe a set of hydro- 
dynamic equations suitable to describe a suspension of 
active particles whose nematic order is allowed to vary in 
space and time as a consequence of activity and flow. 

Let us consider a concentration c of rod-like active par- 
ticles of length £ and mass M suspended in a solvent of 
concentration p so i von t- The total density of the system 
p = Mc + ^solvent is conserved and the fluid is incom- 
pressible. Since the total number of particles is also con- 
stant, the concentration c obeys a continuity equation of 
the form: 

d t c = -V • [c(v + v°) - DVc] , (3) 

where v is the bulk flow velocity, v a is the velocity at 
which the particles actively move relative to the flow, and 
D is the diffusion tensor, which two-dimensional uniaxial 
nematics reads: 

Dij = DoSij + I) (I , , (4) 

where D Q = (D\\ + D ± )/2, D 1 = D|| - D± and D\\ and 
D± are respectively the bare diffusion coefficients along 
the parallel and perpendicular directions of the director 
field. The active current j a = cv a has been modeled 
in different ways and depends on whether the system is 
in a nematic or polar phase. In polar systems, active 
particles are collectively propelled in the direction of the 
macroscopic polarization P; thus v° = t> P with vq the 
average velocity of an individual active particle. Thus, 
for example, for bacterial suspensions vq is constant and 
represents the average swimming velocity of an individ- 
ual bacterium [1], while for mutually propelled particles, 
such as cytoskeletal filaments pushing against each other 
through the action of motor clusters, vq depends on the 
average concentration of the filaments. Thus vq = a\c 
with o.\ = u £ 2 , where Uq is the propulsion velocity for 
unit concentration and is proportional to the rate of ATP 
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FIG. 1: An example of the active currents resulting from Eq. 
(5) in the presence of large distortions of the director field, 
such as those which occur near a disclination. As noted in 
Ref. [10], the tilt in the director field surrounding a +1/2 
disclination (left) results in a collective drift of particles in 
the direction indicated by the red arrows (i.e. toward the 
"nose" of the defect). On the other hand, a —1/2 disclina- 
tion (right) will produce the same amount of incoming and 
outgoing currents and thus zero net flux. 



consumption. This leads to an active current of the form 
j a = aic 2 P [14, 40]. 

For a nematic suspension, on the other hand, the ac- 
tive particles move along n and — n at the same rate; 
thus if the director field is uniform across the system, 
there will be no net flux of particles across an arbitrarily 
small domain and j a = 0. However, in the presence of a 
non-uniform director field or equivalently a non-uniform 
nematic order parameter, there will be regions of fluid 
moving faster than others and thus a current. Such a 
current must depend on the derivatives of the nematic 
tensor rather than on Q^ itself. The simplest term of 
this type with the correct tensorial structure is given by 
vf = voidjQij, which for mutually propelled particles 
gives a current: 



-ax(?djQij , 



(5) 



where ot\ is a constant with dimensions of inverse time. 
The negative sign in Eq. (5) reflects the fact that the 
flux of active particles is directed from regions populated 
by fast moving particles to regions of slow moving parti- 
cles. The active current j° has been derived in the form 
(5) by Ahmadi et al. starting from a microscopic model 
of filaments interacting through a motor cluster [14] and 
later by Lau and Lubensky for swimming bacteria [28] 
(the dependence on the concentration c is different in the 
latter case because of the different microscopic model). 
Ramaswamy and coworkers argued that such active cur- 
rents are responsible for the existence of giant fluctua- 
tions in the number density of active nematic particles 
in the presence of noise [10, 17]. Since spatial variations 
of the order parameter were neglected in that work, the 
only driving force for active currents arose from curva- 
ture (i.e. tilt) in the director field orientation, hence the 
name "curvature- induced currents" coined in [17]. Here 
we show that a more complete description, in which both 
the orientation of the director field and the nematic or- 



der parameter are allowed to vary, leads to additional 
complex phenomena. 

Next we construct a set of hydrodynamic equations 
for the nematic tensor Qij. These can be written in the 
generic form: 



[St 



V]Q 
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(v) 



(a) 



(6) 



where the rates , Clff and embody respectively 
the relaxational dynamics, the coupling with the flow, 
and the active contribution to the dynamics of the ne- 
matic tensor. Following Olmsted and Goldbart [41], the 
rates on the right-hand side of Eq. (6) can be obtained 
phenomenologically by constructing all possible traceless- 
symmctric combinations of the relevant fields of the the- 
ory. These are the strain-rate tensor Uij — ^(diVj+djVi), 
the vorticity tensor Wy = ^{piVj — djVi), and the molec- 
ular tensor fly = —SF/SQij defined from the two- 
dimensional Landau-de Gennes free energy F [39]. In 
two dimensions this reads [55]: 

F = J dA[\AQ t] Q %] + \C (QijQij) 2 + \K diQ jk diQ jk ] 

(7) 

where the coefficients A and C characterize the location 
of a second order phase transition. Since tr(Q 2 ) = S 2 /2, 
at equilibrium one has that S = y/—2A/C. In a hard-rod 
fluids, when the IN transition is driven only by the con- 
centration of the nematogens, one can chose for instance: 
A/K = (c* - c)/2 and C/K = c, so that: 



S=y/1- C*/C 



(8) 



Thus for c> c*, S ~ 1 while for c < c* , S = 0. The 
last term in the free energy expression is the Frank elastic 
energy in the one elastic constant approximation [39] . As 
for passive liquid crystals, the relaxational dynamics of 
Qij is driven by the molecular tensor Fl^: 

^ = 7 - 1 i/ u - - + ±S 2 C)Q tJ - KAQij] (9) 

where 7 is a type of rotational viscosity. The coupling 
between nematic order and flow is found by constructing 
all possible symmetric traceless tensors from the products 
of u«, 



and Q^ [41]. This yields: 

&ij = ftWjj + foiWikQkj — Qik^kj) 



+ft 



UikQkj + QikUkj - 2 tr(uQ) Sij 



(10) 



In two dimensions the last term is identically zero and 



simplifies to: 



^ij = Pl u ij + PziWikQkj — Qik^kj) 



(11) 



where the coefficients /3i and (3% can be found by com- 
paring this expression with the standard Ericksen-Leslie 
theory [39, 42], which gives ft = AS and ft = -1 [41]. 
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Here A is the flow-aligning parameter which dictates how 
the director field rotates in a shear flow. In passive liquid 
crystals, for |A| > 1, the director tends to align to the flow 
direction at an angle 9q such that cos 26*o = 1/A, while for 
|A| < 1, it forms rolls across the system. These regimes 
are known as "flow aligning" and "flow tumbling" re- 
spectively. The value of A has even greater significance 
in active systems; together with the magnitude of forces 
exerted by the active particles it dramatically influences 
the flow behavior and rheological properties of the system 
[19, 29]. 

The active contribution to the dynamics of the nematic 
tensor was derived by Ahamadi et al. [14] and is directly 
proportional to the nematic tensor 

fig } = aoQij . (12) 

However, given the structure of Eqs. (6) and (9), such a 
term can be simply incorporated into the molecular field, 
leading to a redefinition of the critical concentration c*. 
Here we will drop this explicit dependence for sake of 
brevity, but remember that that the critical concentra- 
tion c* associated with the IN transition does depend 
on the activity. The corresponding reactive stress ten- 
sor can be obtained from Eqs. (6), (9) and (11) using 
the standard energy conservation and entropy produc- 
tion argument (see, for example, Ref. [42]). This gives, 
after some algebra: 

= —\SHij + QikHkj — HikQkj (13) 

Finally the flow velocity obeys the Navier-Stokes equa- 
tion, with the total stress tensor given by: 

cry = 2i]Uij - p8 tJ + a\f + a 2 c 2 Q l: j (14) 

where r/ is fluid viscosity of the fluid and p the pressure. 
The last term was established in the seminal work of Pcd- 
ley and Kessler [1] and represents the tensile/contractile 
stress exerted by the active particles in the direction of 
the director field n. The c 2 dependence again arises be- 
cause the propulsion, in our analysis, is provided by pair 
interactions of the filaments through motors. A detailed 
derivation can be found in Ref. [40]. 

Summarizing, the hydrodynamics of an incompress- 
ible nematic suspension of mutually propelled filaments 
is governed by the following set of differential equations 
for the particle concentration c, the nematic tensor Qij, 
and the flow velocity v (with components Uj): 

pd t v l = rjdfvi - d t p + djT i3 (15) 

[d t + Vidi]c = di[(D Sij + DiQij)djC + a±c 2 djQij] 

[d t + Vidi]Qij = XSuij + QikHkj - u lk Q kj + 7 _1 i?ij 

where we defined r^- = + ai(?Qij. Here we have ne- 
glected the inertial convective term in the velocity equa- 
tion because we are interested in fluids of cytoskeletal fila- 
ments and motor proteins for which the typical Reynolds 



number is small. However, since both the nematic order 
and the concentration of particles can be advected by the 
flow, the associated terms in the equations for c and Q irj 
cannot be neglected. 

The dynamics of such an active nematic suspension 
is governed by the interplay between the active forcing, 
whose rate r" 1 is proportional to the activity parameters 
ai and ct2, and the relaxation of the passive structures, 
the solvent and the nematic phase, in which energy is dis- 
sipated or stored. The response of the passive structures, 
as described here, occurs at three different time scales: 
the relaxational time scale of the nematic degrees of free- 
dom f 2 /(7 _1 if), the diffusive time scale £ 2 /D , and the 
dissipation time scale of the solvent pL 2 /rj, where L is the 
system size. While the presence of three dimensionless 
parameters makes for a very rich phenomenology, for sim- 
plicity we choose parameter values in this work so that 
the three passive time scales are of the same magnitude 
r p . When r a 3> r p , the active forcing is irrelevant and the 
system behaves like a traditional passive suspension. On 
the other hand, when r a ~ r p , the passive structures can 
balance the active forcing leading to a stationary regime 
in which active stresses are accommodated via both elas- 
tic distortion and flow. Finally, when r a -C t p the pas- 
sive structures respond too slowly to compensate active 
forces, leading to a dynamical and possibly chaotic in- 
terplay between activity, nematic order and flow. In the 
rest of the paper, we quantify these different regimes. 

For further manipulations it is convenient to make the 
system dimensionless by scaling all lengths using the rod 
length £, scaling time with the relaxation time of the 
director field r p = £ 2 /(7~ 1 ii"), and scaling stresses by 
the elastic stress a = Ki~ 2 . 



III. CHANNEL GEOMETRY 

A. Overview 

The simplest geometry in which to analyze the hy- 
drodynamic equations given in the previous section is 
a two-dimensional channel of infinite length and finite 
width. This geometry has been studied in detail for ac- 
tive nematic and polar suspensions under the assumption 
of constant magnitude of the nematic order parameter 
[12, 19-21, 29]. The most striking feature of active ne- 
matic fluids in a channel is the "spontaneous flow transi- 
tion" : when the activity parameter a 2 is increased past 
a threshold, the system goes from a stationary state in 
which the director field is parallel to the walls of the chan- 
nel to a state of non-uniform orientation and flow. Here 
we show that lifting the assumption of constant nematic 
order parameter leads to a second transition to oscilla- 
tory flow not considered previously in theories of active 
nematics. 

We consider a channel of infinite length along the x 
direction of a Cartesian frame and finite width L along 
y. The channel is bounded by no-slip surfaces at y = 
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and y = L. Assuming translational invariance in the x 
direction, the flow field is completely defined by the ve- 
locity field v x — v x (y),v y — since the incompressibility 
condition implies that 

V • v = dyv y = 0; (16) 

thus, v y = since the fluid is confined. The strain rate 
tensor has only one non-zero component u xy = d y v x /2 
1 1 2, Calling 9 the angle between the director field n and 
the x axis, the nematic tensor can be expressed in the 
simple form: 

5 / cos26» sin 29 \ , . 

Q ~ 2 I, sin20 -cos26> ) ^ 

Under these conditions the hydrodynamic equations (15) 
simplify to: 

pd t v x = d y cr X y 

d t c = d y [(D - \D l S 2 cos29) d y c- \aic 2 ] d y (S cos 26*)] 

d t S = d^S - S[A{c) + \C{c)S 2 - Xu sin 29 + 4(d y 9) 2 } 

d t 9 = d 2 y 9 + 2S- 1 d y S d y 9 - ±u(l - \ cos29) (18) 

where we have assumed that the inertia of the active par- 
ticles is negligible and the dependence of the coefficients 
A and C on the concentration c is as explained previously 
(The procedure for decoupling the angle 9 and the order 
parameter S is given in Appendix A). Furthermore, the 
total shear stress is finally given by: 

a xy = rju+ \a 2 c 2 S s\n29 + ±S\shi29d 2 S 

+25(1 - A cos 29)d y S d y 9 + 25(1 - A cos 20)9*0 

+SX sin 29[A{c) + \C{c)S 2 + A{d y 9) 2 ] . (19) 

To complete the formulation of the problem, we need to 
specify some boundary conditions. Here we assume that 

v x (0) = v x (L) = 0, 0(0) = 0(£) = 0, 

c'(0) = c'(L) = 0, S'(0) = S'(L) = 0. (20) 

Here the condition on 9 assumes strong anchoring with 
the filaments parallel to the wall at the boundaries, and 
the conditions on c and S imply that there is no current 
flowing through the walls. In particular, from Eq. (3) 
and (5): 

j y = (D - iDxS 2 cos 2(9) d v c - \a lC 2 d v {Scos 29) . 

(21) 

The condition j y (0) — j y (L) — requires d = and S' — 
at the boundaries. As initial conditions we take c = c 
(constant), v x = and 9 and S randomly distributed. 

We solved the hydrodynamical equations (18-19) with 
the boundary conditions (20) numerically after dropping 
all derivatives of order higher than two in the equation 
for v x to avoid the use of fictitious boundary conditions. 
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FIG. 2: The hydrodynamic fields c, S, 9 and function 
of y/L in the spontaneously flowing regime for the channel 
geometry obtained by solving Eqs. (18) with boundary con- 
ditions (20). The initial concentration of the active particles 
is set to co = 2c*, »2 = 1.5 and ai = 0.1ct2- The other 
parameters are A = 0.1, L = 5 and r\ = Do = D\ = 1. 
For these values, the spontaneous flow transition occurs at 
a 2 = 1.24069. 

In all our numerical calculations we set ot\ — 0.1a2, al- 
though other parameters are changed as indicated. Our 
simulations show that the system exhibits three different 
regimes determined by the values of the activity param- 
eter CV2 and the flow-alignment parameter A. For small 
activity, the homogeneous stationary state is the only 
stable solution, with: 

v x =0, c = c , S=y/l-<*/c, 9 = (22) 

Upon increasing a 2 and taking < A < 1, the system 
undergoes a transition to a steady state in which 9 and 
S vary across the system. In addition the flow velocity v x 
is non-zero and reaches its maximum in the center of the 
channel. Fig 2 shows a plot of the hydrodynamic fields as 
a function of y/L. We note that the spatial variations of 
the order parameter S are not localized, thus this regime 
is equivalent to the spontaneously flowing state identified 
in earlier studies of active nematics [12, 19-21, 29]. It is 
worth noting that, in this steady state, the nematic order 
parameter is anti-correlated with the concentration. This 
feature, which might appear counterintuitive in compar- 
ison to the passive case, is a non-equilibrium effect that 
arises due to the balance between diffusive and active 
currents (Eq. 5). Assuming n uniform, the total particle 
current is given by: 

j y - -DdyC- a lC 2 d y S (23) 

Since stationary solutions in this geometry require j y = 
0, the active current out of regions with large S is bal- 
anced by diffusion from regions with large c. 

Upon increasing a 2 the system undergoes a further 
transition to a regime in which the order parameter S, 
the tilt angle 9, and the velocity oscillate in time, with a 
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FIG. 3: The hydrodynamic fields c, S, 9 and v x in the center 
of the channel as a function of time in the oscillatory regime 
obtained by solving Eqs. (18) with boundary conditions (20) 
for Q2 = 3 and the other parameters as in Fig. 2. 



frequency that increases with a 2 . In Fig. 3 we illustrate 
the oscillatory behavior by showing a plot of the hydro- 
dynamic fields c, S, 9 and v x in the center of the channel 
(y = L/2) as a function of time. Both the spontaneously 
flowing and the oscillatory regime occur in the nematic 
phase, when the concentration c > c* . For c < c* , on 
the other hand, the isotropic homogeneous state with no 
flow is the only solution. 



B. Linear stability analysis 

To understand the result of our numerical simulations 
and the onset of spontaneous flow we turn to stability 
analysis of the base state. Letting ip — {c, 5, 9, v x }, we 
consider: 



¥>(j/>*) = ¥>o + e<Pi(y,t) 



(24) 



with (fo = {co, So, 0, 0} the stationary homogeneous so- 
lution and e <C 1. Substituting this ansatz into the hydro- 
dynamic equations (18) yields a linearized system that 
may be written in block-diagonal form as: 



A 
B 



with: 



A = 



{D -\D x S Q )d 2 y ~\a lC 2 d 2 y 
|5 (l-5 2 ) ~c S 2 + d 2 



and: 



B 



d 2 

v 

a 2 clS Q d y 



i(i -AH 

Vd 2 



(25) 



(26) 



(27) 



a 2 




A 



FIG. 4: Phase diagram of the flow behavior in the channel 
geometry, presented in the (A,«2) plane (with ai = 02/10). 
The boundary line separating the stationary state (5) from 
the steady flow state (SF) is given by Eq. (31). The phase 
boundary of the oscillatory regime (OF) was obtained nu- 
merically. Other parameter values are r\ — Do = D\ — 1, 
cq = 2c*, and L — 5. 



operator. To calculate the critical value of «2 we must 
solve the homogeneous system: 



B n d% 



Bi 2 d v v 1 = 



B 2 2d 2 Vl + BzLdyBi = 



(28) 



with the boundary conditions (20). This implies that the 
only possible forms for 9\ and v\ are: 



n ' l' 27rn 



vi = C 2 



1 — cos 



'27rn 



Substituting these forms into (28) yields 



(27m) 



t(l-A) 



C 2 = 



2ima2C a Sn 
L 



(2im) 2 T] 



c 2 = o 



(29a) 



(29b) 



(30) 



which together with the requirement for C\ and C 2 to be 
non-zero yields the following critical value of a 2 : 

* 8yTT 2 n 2 
c^L z S (l - A) 

We thus see the first unstable mode corresponds to n = 
1, which along with the critical value a 2 is consistent 
with that seen in our numerical simulations. The phase- 
diagram in Fig. 4 summarizes the flow behavior for the 
channel geometry. 



IV. PLANAR GEOMETRY 



A. Overview 



The spontaneous-flow instability is triggered by the cou- 
pling between orientation and flow embodied in the B 



We now turn to the case of an active nematic fluid in a 
two-dimensional square domain with periodic boundary 



conditions. We numerically integrated the hydrodynamic 
equations of Sec. II using a vorticity /stream- function 
finite difference scheme on a collocated grid of lattice 
spacing Ax — Ay — 0.078. The time integration was 
performed via a fourth order Runge-Kutta method with 
time step At — 10~ 3 . As illustrated in Sec. IV B, the 
vorticity /stream- function method requires one to solve a 
Poisson equation at each time step in order to calculate 
the two components of the flow velocity. This was per- 
formed efficiently with a P-cycle multigrid algorithm [43] . 
As initial configurations we considered a homogeneous 
system where the director field is aligned along the x axis 
and subject to a small random perturbation in density 
and orientation. Thus c = c + e, 6 = e, S = \J\ — c* /c 



10 



and v 7 



0, where e is a random number of zero 



mean and variance (e ) = 10 . The equations were 
then integrated from t = to t = 10 3 , corresponding to 
10 6 time steps. Except where mentioned otherwise, the 
numerical calculations described in this section use the 
parameter values a.\ =0,2/2, r\ = Dq = D\ = 1, A = 0.1, 
c = 2c* and L = 10. 

At low activity, the system relaxes quickly to a sta- 
tionary homogeneous nematic state with: 



V X = Vy = , 



CO , 



S=y/1-C*/C 



= 0. 



Upon raising the activity above a critical value a|, with 
af « 0.37 for the parameters of our calculation, this state 
becomes unstable to a flowing state. The behavior of the 
spontaneously flowing solution, in this two-dimensional 
periodic domain, is substantially different than the quasi- 
one-dimcnsional system discussed in Sec. III. For values 
of Cf2 slightly above a\ ~ the system divides into two 
bands flowing in opposite directions. The direction of the 
streamlines is dictated by the initial conditions which, in 
this case, favor a flow in the x direction. Moreover the 
solution is constant along the flow direction (see Fig. 5). 

The structure of the bands can be inferred from the 
plots in Fig. 6 showing the various hydrodynamic fields 
along the y direction. The yellow region indicates the 
extent of a band. Both the flow velocity and the concen- 
tration are maximal at the center of a band. The max- 
imum in the velocity, in particular, is associated with 
a very sharp variation in the orientation of the director 
field (see the bottom- left panel of Fig. 6). This rapid 
variation of the director field generates a large elastic 
stress, which is balanced by the release of viscous stress 
through the increase in the local flow velocity. The ne- 
matic order parameter, on the other hand, is minimal in 
the center of a band due to the balance between diffusive 
and active currents discussed in Sec. III. As in the case 
of spontaneous flow in the channel geometry, here too the 
variations in concentration and the order parameter are 
relatively small and not localized. 
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FIG. 5: The velocity field (top) and the director filed (bot- 
tom) are superimposed on density plots of the concentration 
(top) and the nematic order parameter (bottom) for «2 = 0.4 
obtained by solving Eqs. (15) with periodic boundary condi- 
tions. The colors indicate regions of large (green) and small 
(red) density and large (blue) and small (brown) nematic or- 
der parameter. The flow consists of two bands traveling in 
opposite directions. The director field is nearly uniform in- 
side each band. Parameter values are Co = 2c*, ai — 02/2, 
A = 0.1, L — 10 and 77 = Do = Di — 1. For these values, the 
spontaneous flow transition occurs at af = 0.37 



B. Linear stability analysis 

To understand these behaviors, we analyze the linear 
stability of the stationary homogenous state in the two- 
dimensional periodic domain. In order to ensure the in- 
comprcssibility condition V • v = it is convenient to 
rewrite the Navier-Stokes equation in terms of vorticity 
and stream function, by writing: 



d y ip v y = -d x ip 



(33) 



so that the incompressibility condition is automatically 
satisfied and the vorticity field is given by: 



2iO X y 



8 X Vy 



8yV X 



(34) 




FIG. 6: Hydrodynamic fields c, S, 9 and v x as a function 
of y obtained by solving Eqs. (15) with periodic boundary 
conditions. The parameter values are the same as in Fig. 5. 
The yellow region indicates the extent of a band shown in 
Fig. 5. 



The two-dimensional Navier-Stokes equation can be ex- 
pressed in terms of u; by: 



d t uj = r]Aw + d 2 x T yx + d xy r yy - d yx r xx 
where we defined: 



y'xy 



Tij — —XSHij + QikHkj — HikQkj + OL2(?Qi 



(35) 



(36) 



From Eq. (34) we see that the stream-function ip is re- 
lated to the vorticity lu through a Poisson equation of the 
form: Atp = —u. Consistent with the numerical calcu- 
lations, we consider a nearly uniform suspension of ne- 
matogens whose director field is approximatively aligned 
along the x direction. Thus c(x, t) — cq + eci(x,i) and 
n(x, t) = x + e n.i(x, t). Analogously, the nematic tensor 
can be expressed to first order in e as: 

S 

ftj( x , *) = irikxSjx ~ SiySjy) + e Q\f (37) 



with Sq = yl — c* Jcq . As in the quasi-one-dimensional 
case, we use the compact notation cp = {c, Q xx , Q X y, lo} 
and write the perturbative expression: 



p(x,t) = ¥> (0) +e¥J (1) (x,t) 



(l) 



(38) 



To enforce periodic boundary conditions on a square do- 
main, we look for solutions of the form: 



¥> (1) (x,*) 



oo oo 
n— — oo m=—oo 



(nx-\-my) 



(39) 



The Fourier components of the stream-function are re- 
lated to those of the vorticity by: 



(40) 



With this choice the linearized hydrodynamic equations 
reduce to a set of coupled of linear ordinary differential 
equations for the Fourier modes tp nm : 



(41) 



with the matrix A nm given in Appendix B. The first 
mode to become unstable is the transverse excitation 
(n, m) = (0, 1) associated with the block-diagonal ma- 
trix: 



aoi 
doi 



with: 



aoi 



-^(2D -£»i5 ) 



L 2 



4c 



S 



4tT_ 
L 2 



(42) 



(43) 



(44) 



The instability first arises from the coupling between lo- 
cal orientations and flow (unless ot\ 02) ■ The critical 
value of «2 is obtained by examining the eigenvalues of 
the matrix doi given by: 




A., 



2tt 2 (1 



V- 

-4n 2 S 2 (l 



77) ± 7T r 



L 2 
A) 2 



2n 2 {l-r 1 ) 2 
a 2 clS L 2 (l 



A) (45) 



When the real part of the above eigenvalues becomes 
positive, an instability ensues: this corresponds to a 2 
larger than the critical value: 



^ 2 [2 V + S 2 (l - \f 
c 2 L 2 S (l-X) 



(46) 



The origin of the instability of the homogeneous station- 
ary state is the same for the ID channel and the 2D 
domain and is related to the interplay between the orien- 
tation of the director field and the shear flow driven by 
the internal active stresses. To illustrate this point let 
us consider a two-dimensional nematic fluid in a station- 
ary state with the director field aligned, say, along the 
x axis of an arbitrary reference frame. The active stress 
produced by the action of the motors powers a collective 
motion of the nematogens. However, the director field 
rotates in the presence of shear flow for A 7^ 1, which 
generates elastic stress. For small activities, the elas- 
tic stiffness dominates and suppresses flow, while above 
the critical value of a| activity dominates and drives col- 
lective motion. Higher levels of nematic order focus the 
sources of active stress and thus require lower activity lev- 
els to destabilize the homogenous stationary state (lower 
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FIG. 7: Hydro dynamic fields c, S 9 and uj at the center of 
the box as a function of time obtained by solving Eqs. (f 5) 
with periodic boundary conditions for a2 = 1.5 and the other 
parameters as in Fig. 5. 



In a "dry" system (i.e. v x — v y — and «2 = 0) the 
homogeneous state becomes unstable solely as a conse- 
quence of the coupling between density and orientation 
fluctuations expressed by the matrix a nm . In this case, 
the first modes to become unstable are the transverse 
mode (1,0) and the longitudinal mode (1,0) associated 
with the matrix (see Appendix A): 



OlO 



-fr (2D + D 1 S ) 
' So 



4c 



4n A 
L 2 

-c S 2 



£*iCo 

4tt 2 
L 2 



(47) 



Simple algebraic manipulations can be used to show that 
the real part of the eigenvalues of aoi and a%o becomes 
positive when a± is larger in magnitude than the critical 
value: 



a-, 



2(±2D -D 1 S Q ){c L 2 S$+4:W- 
c c*S a L 2 



(48) 



where the plus sign correspond to the (1, 0) mode and 
the minus to the (0, 1) mode. This instability, which oc- 
curs in absence of hydrodynamics, has been described in 
various contexts (see for example [22, 23] and references 
therein). We refer the reader to these works for a de- 
tailed discussion while in the rest of this article we focus 
on hydrodynamic phenomena. A thorough discussion on 
the instability of the homogeneous state in "dry" and 
hydrodynamic systems can be found in [25] . 



C. Relaxation Oscillations 

Upon increasing the activity parameter a 2 above a sec- 
ond critical value (with a\ « 0.41 for our default pa- 
rameter values), the spontaneously flowing state evolves 
into a pulsatile spatial relaxation oscillator. Fig. 7 shows 
a plot of the various hydrodynamic fields as a function of 
time for oti = 1.5. In this regime the dynamics consists of 
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FIG. 8: Dynamics of an active "burst" for the trajectory 
shown in Fig. 7, with a.2 = 1.5. The flow velocity at the 
point x = y = L/3 is shown as a function of time over the 
course of a director field rotation (top left) and the director 
field is shown for the three labeled time points. Between two 
consecutive bursts the system is homogeneous and uniformly 
aligned. During a burst, nematic order is drastically reduced 
in the whole system and the director undergoes a distortion 
with a consequent formation of two bands flowing in opposite 
directions. After a burst, a stationary state is restored with 
the director field rotated of 90° with respect to its previous 
orientation. 



a sequence of almost stationary passive periods separated 
by active "bursts" in which the director switches abruptly 
between two orthogonal orientations. During passive pe- 
riods, the particle concentration and the nematic order 
parameter are nearly uniform across the system, there is 
virtually no flow, and the director field is either paral- 
lel or perpendicular to the x direction. Eventually this 
configuration breaks down and the director field rotates 
by 90° (see Fig. 8). The rotation of the director field 
is initially localized along lines, generating flowing bands 
similar to those discussed in Sec. IV A. The temporary 
distortion of the director field as well as the formation of 
the bands is accompanied by the onset of flow along the 
longitudinal direction of the bands. The flow terminates 
after the director field rotates and a uniform orientation 
is restored. The process then repeats. 

Remarkably, the rotation of the director fields occurs 
through a temporary "melting" of the nematic phase. As 
shown in Fig. 7, during each passive period the nematic 
or der para meter is equal to its equilibrium value So = 
y/l — c* Jc (S = l/v2 because of the choice of cq = 2c*), 
but drops to ~ giSo during rotation. The reduction of 
order is system-wide, but, as shown in the bottom-left 
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FIG. 9: (Left) The vorticity u, the orientation of the director 
cos#, and the nematic order parameter 1 — S/So are shown 
for the point x — y = L/3 over the course of a burst for the 
trajectory shown in Fig. 7 with ai = 1.5. The data is from 
the numerical integration and the vorticity is normalized so 
that its maximum value is one. (Right) A close-up of the 
same data during the onset of a burst. 
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FIG. 10: (Left) The average nematic order parameter (S) = 
J dA/L 2 S(x) and the total shear stress a xy are shown over 
several bursts for for the trajectory shown in Fig. 7 with 02 = 
1.5. (Right) The frequency of bursts is shown as a function 
of Q2 with other parameters as in Fig. 7. 



panel of Fig. 8, is most pronounced along the boundaries 
between bands. Without this transient melting (i.e. if 
the magnitude of S is not allowed to vary) , the distortions 
of the director field required for a burst are unfavorable 
for any level of activity. 

A closer look at the dynamics of an individual oscilla- 
tion elucidates the mechanism of the instability. Fig. 9 
shows the flow (represented by the vorticity), orientation, 
and the nematic order parameter as a function of time 
for o?2 = 1.5. Beginning from the homogeneous state, 
the active forcing generates a gradual increase in flow, 
and the system evolves in a manner similar to that of 
the spontaneous flow regime described in section IV A. 
As described there, the resulting shear flow causes the 
nematic director to rotate, generating elastic stress that 
competes with the active stress. Above the critical value 
of cv2, however, the elastic stress is never sufficient to bal- 
ance the active stress and the banded flow configuration 
becomes unstable to melting of the nematic phase. Im- 
portantly, the instability occurs only once the flow and 
director rotation have reached a threshold level; thus, 
there is a significant delay during which the nematic or- 
der parameter is nearly constant. Once melting occurs, 
the stress is rapidly released during reorientation. The 
timescale of the oscillation is given by the time required 
for the flow and director rotation to reach their threshold 
values, and thus decreases with an increase in a 2 above 

The physical origin of the oscillatory dynamics in our 
model of active nematic suspension has to be ascribed to 
the existence of multiple time scales in a system that is 
internally driven. As we mentioned in Sec. II, one time 
scale is set by the rate at which the active forcing oc- 
curs and is r a = 77/(0:200). A second time scale is related 
to the relaxational dynamics of the fluid microstructures 
(i.e. the solvent flow field, the director field, and the ne- 
matic order parameter) and is given by t p = ^/(■y~ 1 K) 
(the time unit in all numerical and analytical calcula- 
tions). When the two time scales are comparable, the 
active forcing is accommodated by the microstructures 
leading to a distortion of the director field and a steady 
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FIG. 11: (Left) A typical trajectory of the variable Q from 
Eqs. (51) for a slightly above the critical value |??(2a + rjk 2 ). 
(Right) The same limit cycle in the (Q,«)-plane. The black 
dashed line is the ii = nullcline and the black solid line is 
the Q — nullcline. 



flow. However, when the active forcing occurs at a larger 
rate the microstructures fail to keep up, revealed above 
by the instability to melting of the nematic phase. This 
lag results in oscillatory dynamics and eventually chaos. 
Similar oscillatory phenomena have been found in mod- 
els of complex fluids under shear. Cates and cowork- 
ers discussed specifically the effect of a slow response of 
the microstructure to an external shear and showed how 
such a phenomenon can be naturally described via the 
FitzHugh-Nagumo equation [44-46]. 

To illuminate the origins of the relaxation oscillations 
we construct a simplified version of the hydrodynamic 
equations that retains the minimal features required to 
exhibit oscillatory phenomena: the coupling between ac- 
tive forcing and the fluid microstructure and the variable 
nematic order. The purpose of the following calculation 
is not to rigorously analyze Eqs. 15, but rather to iden- 
tify basic physical mechanisms that can drive oscillations 
and to illustrate the effect of different timescales in the 
system. 

Let us then consider the following simplified version 
of the hydrodynamic equations for the quantities Q xy 
and u xy which represent respectively the liquid crystal 
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degrees of freedom and the flow field. 

Qxy V"xy 7 H- X y , (49a) 
ii xy = A(r)u xy + aQ xy ) , (49b) 

obtained by treating c and Qxa; as constants and by sim- 
plifying the coupling between the nematic tensor and 
flow, as compared to to the complete phenomenological 
construction discussed in Sec. II. Here, variations in the 
nematic order parameter are embedded in the Landau-de 
Gennes free energy within H xy . Moving to Fourier space, 
Eqs. (49) can be rearranged in the form: 

Qxy = u xy + 7 - 1 [(l^l - CQL - k 2 )Q xy - CQly] , 

ii xy = -k 2 (r\u xy + aQ xy ) , (50) 

Finally, by taking Q = Q xy , u = -u xy , a = 7 _1 (|^4| - 
CQ XX — k 2 ) and b = 7~ 1 C, one obtains: 

Q = aQ-bQ 3 -u (51a) 

u = k 2 (aQ — rju) , (51b) 

equivalent to the spatially homogeneous FitzHugh- 
Nagumo model or the generalized van der Pol oscil- 
lator [47, 48]. In our periodic square domain k 2 = 
(2nir/L) 2 + {2mn / ' L) 2 , with n and m integer numbers. 
The nullclines of the dynamic system (51) are given by: 

u = Q(a - bQ 2 ), u = aQ/rj (52) 

There are, in general, three fixed points P = (Q, u): 

(53) 

For a < r](2a + r]k 2 )/3 the origin Pq is a saddle point, 
while P± are stable nodes. A trajectory starting from an 
arbitrary [Q 1 u) point will then converge to a stable state 
characterized by a finite strain-rate that matches the ac- 
tive stress aQ: rju = aQ = a^J (a — ajrjfjb. This fixed 
point represents the usual spontaneously flowing state 
(Fig. 5). For a > r](2a + r/fc 2 ), P± become unstable and 
the system exhibits relaxation oscillations. Fig 11 shows 
a typical trajectory and a phase-plane plot showing the 
flow of trajectories in u and Q space. In this regime the 
dynamics consists of slow relaxations, when the trajec- 
tory is close to the cubic nullcline (Q = 0), interspersed 
with rapid large jumps in Q when the trajectory reaches 
the unstable portion of the cubic nullcline. By inspec- 
tion of Eq. (51b) the frequency of the oscillations is 
given by v ~ k 2 a. To expand on the assertion that relax- 
ation oscillations arise when the passive timescales exceed 
that of the active forcing, the critical active rate can be 
obtained by rewriting a c in terms of the characteristic 
timescales defined above as St^ 1 — (2ar p _1 + ^ 2 fc 2 r d _1 ), 
with r a = rj/a. 



Numerical simulations of the full system (15) exhibits 
a much richer behavior than that captured by Eqs. (51), 
but the qualitative dependence of the dynamics with re- 
spect to «2 persists. The origins of the kink at a-i = 1-35 
are unclear at present, but it does not correspond to ex- 
citation of a spatial mode of larger wave-number. 

It is interesting to study how the three regimes de- 
scribed so far change as the size of the system is increased. 
Fig. 12 shows a phase diagram of the various dynamical 
regimes for the full equations in the plane (L,«2)- Upon 
increasing the size L of the system, the critical value 
of q;2 separating the spontaneous flow and the oscilla- 
tory regime decreases and merges with the lower phase 
boundary [whose expression is given in Eq. (46)] for 
11 < L < 12. Thus we expect that in large samples, 
the instability of the homogeneous state will lead directly 
to oscillatory and then chaotic dynamics. The latter is 
described in the next section. 

It is important to emphasize that the excitability de- 
scribed here for active nematic fluids is a purely hydro- 
dynamic phenomenon that arises as a consequence of the 
existence of multiple time scales in the system, when the 
dynamics of the flow lags with respect to the rate of the 
active forcing exerted at the microscopic scale. This phe- 
nomenon is thus very different from the large scale fluctu- 
ations previously observed in simulations with noise and 
no hydrodynamics [18, 22]. Furthermore, the excitabil- 
ity seen here is quite different from that seen in many 
biological systems where the relaxation oscillations arise 
from heavily regulated networks of chemical and electri- 
cal signals, in contrast with what see in our model where 
they emerge directly from physical interactions among 
the constituent components of an active fluid such as the 
cytoskeleton in a cell. 



D. Chaotic regime 

When the activity a2 is further increased past a third 
critical value a\ 1 with a% ~ 2 for our default parame- 
ters (Fig. 12), the flow becomes chaotic. The route to 
chaos takes place through a disordering of the flip-flop 
dynamics described in the previous section. Initially the 
dynamics is still characterized by periods of low activ- 
ity alternating with bursts during which nematic order is 
temporarily lost and the director field rotates. In Fig. 13 
we show the time course of several hydrodynamic fields 
in a typical trajectory for a 2 = 2.3. 

In this chaotic regime, the structure of the 
flow presents some coherent features typical of two- 
dimensional turbulence. For example, in Fig. 14 we show 
a representative snapshot of the flow velocity superposed 
on the concentration field, and the director field super- 
posed on the nematic order parameter. We see that the 
flow is characterized by large vortices that span the sys- 
tem size, with the director field organized into "grains" of 
uniform orientation separated by grain boundaries that 
span the entire sample. Comparison of the two plots in 



12 




10 



FIG. 12: Phase diagram for the stationary (£), spontaneous 
flow (F), relaxation oscillation (O) and chaotic (CT) regimes in 
the plane (L, 02) for the full equations Eqs. (15) with periodic 
boundary conditions. The dots are obtained from numerical 
integration. The green solid line, separating the stationary 
and flowing state, is given by Eq. (46). The red dashed line, 
separating the spontaneously flowing state and the relaxation 
oscillations regime is interpolated from the numerical data. 
The color gradient at the intersection between the oscillatory 
and the chaotic region, indicates a fuzzy boundary between 
these two regimes. Parameter values are Co = 2c*, ai = 02/2, 
A = 0.1, and rj = D = D 1 = 1 
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FIG. 14: (top) The velocity field superimposed on a density 
plot of the concentration and (bottom) the director field su- 
perimposed on a density plot of the nematic order parameter 
obtained by solving Eqs. (15) with periodic boundary condi- 
tions for 02 = 3 and other parameters as in Fig. 5. The colors 
indicate regions of large (green) and small (red) concentration 
and large (blue) and small (brown) nematic order parameter. 
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FIG. 13: Hydrodynamic fields c, S, 9 and to at the center of 
the box as a function of time obtained by solving Eqs. (15) 
with periodic boundary conditions for Q2 = 2.3 and the other 
parameters as in Fig. 5. 



Fig. 13 reveals that the grain boundaries are the fastest 
flowing regions in the system. Thus the dynamics in 
this regime is characterized by grains with approxima- 
tively uniform orientation that swirl around each other 
and continuously merge and reform, giving rise to a flow 
that appears turbulent. This is similar to other chaotic 
flows in active fluids that have been reported in models 
of dilute bacterial suspensions but which do not include 
liquid crystalline elasticity [49, 50] (also see Ref. [20] for 
a related steady state analysis). 

Fig. 15 shows the energy and enstrophy power spec- 
tra, with the spectral densities E(k) and £l(k) defined 



so that \{v 2 ) = J^dkE(k) and \{lo 2 ) = /" dfcfi(fc) 
are the mean kinetic energy and enstrophy per unit area. 
Although our simple numerical simulations do not span 
a sufficient range of scales to establish any scaling laws 
that are expected of two-dimensional turbulence, there 
are qualitative signatures of such behavior in our model 
of active nematic fluids. We recall that for passive two- 
dimensional fluids, the classic Kraichnan theory of two- 
dimensional turbulence in viscous fluids [51, 52] predicts 
a double cascade through which energy is transfered from 
small to large scales while enstrophy flows from large to 
small scales. At length scales smaller than the injec- 
tion scale, the enstrophy cascade dominates, giving rise 
to energy and enstrophy spectra decaying like fc -3 and 
respectively (modulo logarithmic corrections). The 
fundamental difference between simple viscous fluids and 
the active fluid discussed here is that the forcing acts on 
a molecular scale here, in contrast with the situation in 
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FIG. 15: Energy (top) and enstrophy (bottom) spectra for 
system of size L — 20 obtained by solving Eqs. (15) for Q2 = 
2 and other parameters as in Fig. 5. Dashed lines show 
the graph of the power laws fc~ 3 and A; -1 expected in two- 
dimensional turbulence. 



time scales associated with solvent and liquid crystalline 
degrees of freedom. While we have specifically chosen 
parameter values so that the solvent and liquid crystalline 
degrees of freedom have the same intrinsic timescales, it 
would be interesting to continue the analysis to cases 
with multiple relaxation time scales. 

More generally, the richness of behaviors emerging in 
the present theoretical study of active fluids with liquid 
crystalline order raises an important question: are these 
phenomena observed in real active systems ? And if so, 
how well can hydrodynamic models capture the complex- 
ity of those systems ? In a recent publication Schaller et 
al. [53, 54] reported the observation of many examples 
of the collective dynamics in a motility assay consisting 
of highly concentrated active polar filaments propelled 
by immobilized molecular motors in a planar geometry. 
These include the onset of traveling density bands, oscil- 
latory dynamics in which the average orientation of the 
filaments switches periodically in time, and large scale 
swirling motions. Our results suggest that spatially in- 
homogeneous nematic order is sufficient to drive both an 
oscillatory dynamics of the director field and a swirling 
motion even in the absence of polar order. With this 
work, we hope to have provided a number of testable 
predictions that can be used in combination with exper- 
iments to shed light on the basic physical mechanisms 
governing the dynamics of living or otherwise active mat- 
ter. 



viscous fluids which is forced at scale of the system. This 
suggests that a possible mechanism for turbulence in the 
active fluid described here could involve an inverse en- 
strophy cascade in which vorticity is injected into the sys- 
tem at small scales through the active forcing and then 
transfered to the scales of order the system size. The 
dashed lines in Fig. 15 show the power laws E(k) oc fc~ 3 
and O(fc) cx fc _1 expected for two-dimensional turbulence 
in viscous fluids, which while suggestive are not definitive 
as our numerical methods are inadequate to stringently 
test these ideas quantitatively. However, we hope that 
our simple discussion might serve as a starting point for 
identifying and characterizing active turbulence. 



V. CONCLUSIONS AND OUTLOOK 

In this article we have analyzed in some detail the hy- 
drodynamics of active nematic suspensions in quasi-one 
and two dimensions. By allowing spatial and temporal 
fluctuations in the nematic order parameter, we observed 
a rich interplay between order, activity and flow. Signif- 
icantly, we find that allowing fluctuations in the magni- 
tude of the order parameter S qualitatively changes the 
flow behavior as compared to systems in which S is con- 
strained to be uniform. 

At a minimal level, the behavior of the system can be 
qualitatively understood by comparing the timescale of 
energy input due to activity and the relevant relaxation 
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Appendix A: Nematodynamics via Pauli Matrices 

For some practical application, such as the channel ge- 
ometry described in Sec. Ill, it is desirable to have sep- 
arate hydrodynamic equations for the variables 9 and S, 
rather than having them entangled in the equation for the 
nematic tensor Qy. In two dimensions, this operation 
can be performed rather elegantly by using Pauli matri- 
ces. To see this let us start from the two-dimensional 
nematic tensor expressed in matrix form: 

n _ S f cos 29 sin 29 \ . . 

Q ~ 2 \ sin26» -cos26» ) ^ > 

In order to decouple S and 9, we can introduce the fol- 
lowing matrices: 

cr p = sin 29 <j\ + cos 29 er 3 (A2a) 
7r = cos 29 a i - sin 29 er 3 (A2b) 
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where cr\ and <r 3 are Pauli matrices: 



1 

1 



er 2 



-i 

1 



<7\3 



1 

-1 



The matrices <r p and tt enjoy a number of properties. 
Namely: 



(Tp (Tp — 7T 7T — (5 . 



(TpTV = MT 2 



(A3) 



where 5 is the 2x2 identity matrix. Since the Pauli 
matrices are traceless and Hermetian, so are <r p and 7r. 
An equation for S can be derived straightforwardly by 
expressing: 



n S 



thus: 



dQ 1 fdS' 
2 { ,11 ' 



dt 



Multiplying this expression from the left by <x p and tak- 
ing the trace gives: 

dQ\ 1 fdS\ .„ /d0\ , , dS 



Analogously we have that: 



tr n 



2S 



(A6) 



from which the hydrodynamic equations for S and 9 are 
found in the form: 



dS ( dQ\ 

d0 _ 1 / dQ\ 
d< ~ 25 tr V 



(A7a) 
(A7b) 



Thus, the general hydrodynamic equations of Sec. II can 
be finally recast as follows: 



d t v = V • <r (A8) 
[d t + v-V]S= [\ s u + Qu - uQ + 7 - 1 Jf] £rp 

[d t + v-V}8=^- [Xsu + Qu-uQ + j^H]^ 
[d t + v • V]c = V • [(A)S + £»iQ)Vc + aic 2 V • Q] 



(A4) where we used the notation: [A] a = ti[otA] 



(A5) 



Appendix B: Linearized System 



In Sec. IV B we discussed the linear stability of the 
homogeneous state and we gave an expression for the 
matrix Aoi of the linearized dynamics associated with 
the first unstable mode. Here we give an expression for 
the generic A nm matrix. This can be written in the block 
form: 



(Bl) 



with: 



-^[2D Q (n 2 + m 2 ) + D^n 2 - m 2 )] ^a lC 2 (m 2 

(B2) 

f^So -%(n 2 +m 2 )-c S 2 " 



-^E^aicg 

Km= ( | (B3) 

^^ASn 



(B4) 

^f^(4a 2C 2-Ac*So) ^[a2C 2 ~X(c*~c )S + ^^(n 2 ' ■ ' 
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■ J L V ' 2(n 2 +m 2 ) , 

d nm = | j (B5) 

^a 2 c 2 (m 2 - n 2 ) - ±ff^ (n 2 + m 2 )[n 2 {\ + A) + m 2 (l - A)] -f^(n 2 + m 2 ) 
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